Probing macroscopic quantum states with a sub-Heisenberg accuracy 
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Significant achievements in high-sensitivity measurements will soon allow us to probe quantum 
behaviors of macroscopic mechanical oscillators. In a recent work [arXiv: 0903.0079], we formulated 
a general framework for treating preparation of Gaussian quantum states of macroscopic oscillators 
through linear position measurements. To outline a complete procedure for testing macroscopic 
quantum mechanics, here we consider a subsequent verification stage which probes the prepared 
macroscopic quantum state and verifies the quantum dynamics. By adopting an optimal time- 
dependent homodyne detection in which the phase of the local oscillator varies in time, the condi- 
tional quantum state can be characterized below the Heisenberg limit, thereby achieving a quantum 
tomography. In the limiting case of no readout loss, such a scheme evades measurement-induced back 
action, which is identical to the variational-type measurement scheme invented by Vyatchanin et al. 
but in the context for detecting gravitational waves. To motivate macroscopic quantum mechan- 
ics experiments with future gravitational-wave detectors, we mostly focus on the parameter regime 
where the characteristic measurement frequency is much higher than the oscillator frequency and 
the classical noises are Markovian, which captures the main feature of a broadband gravitational- 
wave detector. In addition, we discuss verifications of Einstein-Podolsky-Rosen-type entanglement 
between macroscopic test masses in future gravitational-wave detectors, which enables us to test 
one particular version of gravity decoherence conjectured by Diosi and Penrose. 



I. INTRODUCTION 

Due to recent significant advancements in fabricating 
low-loss optical, electrical and mechanical devices, we will 
soon be able to probe behaviors of macroscopic mechani- 
cal oscillators in the quantum regime. This will not only 
shed light on quantum-limited measurements of various 
physical quantities, such as a weak force, but also help us 
to achieve a better understanding of quantum mechanics 
on macroscopic scales. 

As a premise of investigating macroscopic quantum 
mechanics (MQM), the mechanical oscillator should be 
prepared close to be in a pure quantum state. To achieve 
this, there are mainly three approaches raised in the lit- 
erature: (i) The first and the most transparent approach 
is to cool down the oscillator by coupling it to an ad- 
ditional heat bath that has a temperature T a( jd much 
lower than that of the environment To. As a result, the 
oscillator will achieve an effective temperature given by 
T e fi = {T j m + T add r add )/(7 m + r add ) with j m and r add 
denoting the damping due to coupling to the environ- 
ment and the additional heat bath, respectively. In the 
strong-damping regime with T add j m , we achieve the 
desired outcome with T e g w T add . Since the typical op- 
tical frequency uiq can be much higher than fcriTo/h, a 
coherent optical field can be effectively served as a zero- 
temperature heat bath. Indeed, by coupling an oscillator 
parametrically to an optical cavity, many state-of-the- 
art experiments have demonstrated significant cooling of 
the oscillator, achieving a very low thermal occupation 
number [l|-|l6|. Similar mechanism also applies to the 
electromechanical system as demonstrated in the exper- 



iments [17H19l |; (ii) The second approach is to introduce 
additional damping via feedback, i.e., the so-called cold- 
damping. The feedback loop modifies the dynamics of 
the oscillator in a way similar to the previous cooling 
case. Such an approach has also been realized experi- 
mentally [20M22| . If the intrinsic mechanical and elet- 
rical/optical qualities of the coupled system arc high, 
those cooling and cold-damping experiments can even- 
tually achieve the quantum ground state of a mechani- 
cal oscillator |23W28l |; (iii) The third approach is to con- 
struct a conditional quantum state of the mechanical os- 
cillator via continuous position measurements. Quantum 
mechanically, if the oscillator position is being contin- 
uously monitored, a certain classical trajectory in the 
phase space can be mapped out, and the oscillator is 
projected into a posteriori state 12911 which is also called 
a conditional quantum state [30|-|35| . Given an ideal con- 
tinuous measurement without loss, the resulting condi- 
tional quantum state of the oscillator is a pure state. 

Recently, we theoretically investigated the third ap- 
proach f or g eneral linear position measurements in great 
details |35(. The analysis of this work is independent 
of the scale and mass of the oscillator — these parame- 
ters will only modify the structure of arising noises. In 
particular, we applied our formalism to discuss MQM 
experiments with macroscopic test masses in future 
gravitational- wave (GW) detectors. We demonstrated 
explicitly that given the noise budget for the design sen- 
sitivity, next-generation GW detectors such as Advanced 
LIGO [36| and Cryogenic Laser Interferometer Observa- 
tory (CLIO) [37} can prepare nearly pure Gaussian quan- 
tum states and create Einstein-Podolsky-Rosen-type en- 
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FIG. 1: (Color online) A schematic plot of a Wigner function 
W(x,p) (left) and the corresponding uncertainty ellipse for 
the covariance matrix V cond (can be viewed as a projection 
of the Wigner function). The center of the plot is given by 
the conditional mean (a; cond ,p cond ). The Heisenberg limit is 
shown in a unit circle with radius given by the zero-point fluc- 
tuation h/(2muj m ). For a pure Gaussian conditional quantum 
state, the area of the ellipse, i.e., 7r det V cond /(2r?ict; m ) 2 , is also 
equal to that of Heisenberg limit. Therefore, the uncertainty 
product detV cond can be served as an appropriate figure of 
merit for quantifying purity of a quantum state. 



tanglement between macroscopic test masses. Besides, 
we showed that the free-mass Standard Quantum Limit 
(SQL) [38l - 440| for the detection sensitivity 
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where m is mass of the probing test mass and f2 is the 
detection frequency, also serves as a benchmark for MQM 
experiments with GW detectors. 

More concretely, a Gaussian conditional quantum state 
is fully described by its Wigner function as shown 
schematically in Fig. [TJ It is given by 



W(x,p) 



1 



: exp 
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2 7 rV / detV cond 

(2) 

Here X = [x - x cond ,p- p coad ] with x cond and p cond de- 
noting conditional means of oscillator position x and mo- 
mentum p, and V cond is the covariance matrix between 
position and momentum. Purity of the conditional quan- 
tum state can be quantified by the uncertainty product 
which is defined as 



U 



■Vdet V c 



TACOI1' 
V XX 



dT^cond 
v pp 



T/'cond'' 
xp 



(3) 



which is also proportional to square root of the area of 
the uncertainty ellipse as shown in Fig. [TJ In Ref. [35j |. 
we related this uncertainty product U of the conditional 
quautum state of test masses in GW detectors to the 
SQL-beating ratio of the classical noise, and the amount 
of entanglement between test masses to the size of the 
frequency window (ratio between upper and lower ends 
of that frequency window) in which the classical noise 
goes below the SQL. 

A state-preparation stage alone does not provide a 
complete test of MQM. This is because the measurement 
data in the state-preparation process only allow us to 



FIG. 2: (Color online) A schematic plot of random walk of 
the conditional quantum state, i.e., its Wigner function, in 
the phase space. Its center is given by the conditional mean 
[a; cond (t), p cond (t)], with the uncertainty given by conditional 
variances V^p PyXp . To verify the prepared conditional quan- 
tum state, the only knowledge that verifier needs to know is 
classical information of the conditional mean provided by the 
preparer if noises are Markovian. 



measure a classical trajectory of the oscillator - quantum 
fluctuations are only inferred from the noise budget, but 
not directly visible. Therefore, the resulting conditional 
quantum state critically relies on the noise model of the 
measurement device. If such noise model is imprecise, it 
will yield severe discrepancies between the actual quan- 
tum state and the conditional one. Therefore, it calls 
for a second measurement stage which has to follow up 
the preparation stage. In this paper, we will address the 
above issue by considering a subsequent state- verification 
procedure, in which we make a tomography of the con- 
ditional quantum state obtained during the preparation 
stage. On the one hand, this verification stage can serve 
as a check of the specific noise model to verify the pre- 
pared quantum state. On the other hand, if we insert an 
evolution stage with the oscillator evolving freely before 
the verification, the quantum dynamics of the oscillator 
can also be probed, which allows us to study different de- 
coherence effects and also check whether a macroscopic 
mechanical oscillator does evolve in the same way as a 
quantum harmonic oscillator or not. 

Since the conditional quantum state undergoes a ran- 
dom walk in the phase space as shown schematically in 
Fig. [2] classical information of the conditional mean, ob- 
tained by the preparer from the measurement data, needs 
to be passed onto the verifier who will then proceed with 
a tomography process. Suppose the state preparation 
stage ends at t = and the preparer obtain a conditional 
quantum whose Wigner function is W(x(0),p(0)). The 
task of the verifier is trying to reconstruct this Wigner 
function by synthesizing marginal distributions of differ- 
ent mechanical quadratures X^(0) from ensemble mea- 
surements at t > 0, and 



A c (0) = i(0)cosC 



P(0) 
muj m 



sinC, 



(4) 



where x(0) and p(0) denote oscillator position and mo- 
mentum at t = and u rn is the oscillation frequency. 
This process is similar to the optical quantum tomog- 
raphy where different optical quadratures are measured 
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FIG. 3: (Color online) A schematic plot of the uncertainty 
ellipses of reconstructed states with the same prepared Gaus- 
sian quantum state but different levels of verification accu- 
racy, which shows the necessity of a sub-Heisenberg accu- 
racy. The center of the plot is given by the conditional mean 
(i co ° d , p cond ). The shaded regimes correspond to the verifica- 
tion accuracy. The Heisenberg limit is shown by a unit circle. 
The dashed and solid ellipses represent the prepared state and 
the reconstructed states respectively. 



FIG. 4: (Color online) Values of reconstructed Wigner func- 
tions on the p — plane, i.e., W IDCOn (x,p — 0), for a single- 
quantum state, obtained at different levels of verification ac- 
curacy. Solid curve shows the ideal case with no verifica- 
tion error. Dashed and dotted curves correspond to the cases 
with a verification error of 1/4 and 1/2 of the Heisenberg 
limit, respectively. The negative regime (shaded) or the non- 
classicality vanishes as the verification error increases. This 
again manifests the importance of a sub-Heisenberg verifica- 
tion accuracy. 



with homodync detections [41J. However, there is one sig- 
nificant difference mechanical quadratures are not di- 
rectly accessible with linear position measurements which 
measure 



x q {t) = x(0) cosw m t + 



P(0) 
muj m 



smw m t, 



(5) 



rather than X^. To probe mechanical quadratures, we 
propose the use a time-dependent homodyne detection 
with the local-oscillator phase varying in time. Given a 
measurement duration of Tint, we can construct an inte- 
gral estimator, which reads 



X 

with 

r-Ti 



lt pfO) 
dtg(t)x(t) oc x(0) cos(' + -^- L 
o mui m 



sinC' (6) 



cos £' 



/ Tint dt g{t) cos uj m t and sin£' 



J lnt dt g{t) smw m t. Therefore, a mechanical quadrature 
X(;> is probed [cf. Eq. Here git) is some filtering 

function, and it is determined by the time-dependent ho- 
modyne phase and also the way how data at different 
time are combined. 

The ability to measure mechanical quadratures does 
not guarantee success of a verification process. In or- 
der to recover the prepared quantum state, it requires a 
verification accuracy below the Heisenberg limit. Phys- 
ically, the output of the verification process is a sum of 
the mechanical-quadrature signal and some uncorrelated 
Gaussian noise. Mathematically, it is equivalent to ap- 
plying a Gaussian filter onto the original Wigner function 
W(x,p) of the prepared state [43|, and thus the recon- 
structed Wigner function is 



W„ 



i{x,p) 



dx' dp' ip(x 



■x',p- 



p')W(x',p') (7) 



where the Gaussian filter i[>(x,p) is given by 



tp(x,p) 



1 



27iVdet V add 



exp 



2^ 



add 



(8) 



with £ = [x,p] and V add denoting the covariance matrix 
for the added verification noise. If the prepared quan- 
tum state is Gaussian, using the property of Gaussian 
integration, the reconstructed Wigner function reads 



W lecon (x,p) 



1 



exp 



2?rVdct V r 
and the covariance matrix V rocon 



_l|- V recon-l|T 
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V 1 ' 



= V 



cond 



radd 



(9) 



(10) 



In Fig. [3l we show schematically the effects of different 
levels of verification accuracy given the same prepared 
conditional quantum state. A sub-Heisenberg accuracy, 
with an error area smaller than the Heisenberg limit, is 
essential for us to obtain a less distorted understanding 
of the original prepared quantum state. In addition, if 
the prepared quantum state of the mechanical oscillator 
is non-Gaussian p^ - |48| . a sub-Heisenberg accuracy is a 
necessary condition for unveiling the non-classicality of 
the quantum state as shown schematically in Fig. 2] and 
proved rigorously in the Appendix [X] 

Verifications of quantum states below the Heisenberg 
limit also naturally allow us to test whether entanglement 
between two macroscopic test masses in GW detectors 
can indeed be established, as predicted in Ref. [H, Hj| , 
and how long such entangled state can survive. Survival 
of macroscopic entanglement can test one particular ver- 
sion of gravity decoherence conjectured by Diosi (49| and 
Penrose (50| . For an individual object, it is not entirely 
clear the classical superposition of what pointer states 
gravity decoherence will drive it into. For an entangled 
state among multiple objects, even though Gaussian, it 
would naturally have to decay into the one that is not 
entangled, within the gravity decoherence timescale. 

As we will show, in order to achieve a sub-Heisenberg 
accuracy, we need to optimize the local-oscillator phase 
of the time-dependent homodyne detection as well as 
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the weight with which data collected at different time 
will be combined. If there is no readout loss, this opti- 
mization automatically will give a detection scheme that 
evades measurement-induced back action, the same as 
the variational-type measurement scheme proposed by 
Vyatchanin and Matsko |42j for detecting gravitational- 
wave signals with known arrival time. Since in a single 
measurement setup, different quadratures do not com- 
mute with each other, namely 



i h 
mu) m 



sin(C-C'), 



(11) 



one needs multiple setups and each makes ensemble mea- 
surements of one particular quadrature with a sub- 
Heisenberg accuracy - the synthesis of these measure- 
ments yields a quantum tomography. 

As a sequence to Ref. [35[ and to motivate MQM ex- 
periments with future GW detectors, we will also focus 
on the same parameter regime where the characteristic 
measurement frequency is much higher than the oscil- 
lator frequency and the oscillator can be treated as a 
free mass. In addition, we will consider situations where 
the spectra of the classical noise can be modeled as be- 
ing white. Non-Markovianity of noise sources - although 
they certainly arise in actual GW detectors [35[ and will 
be crucial for the success of a real experiment — is a 
rather technical issue. The non-Markovianity will not 
change the results presented here significantly, as we will 
show and address in a separate paper [5l| . 

This paper is organized as follows: in Sec. [Til we will 
formulate the system model mathematically by writing 
down the Heisenberg equations of motion; in Sec. IIII1 
we will provide a timeline for a full MQM experiment 
with preparation, evolution and verification stages, and 
use simple order-of-magnitudc estimates to show that 
this experimental proposal is indeed plausible; in Sec. |V1 
we will evaluate the verification accuracy in the pres- 
ence of Markovian noises (largely confirming the order-of- 
magnitude estimates, but with precise numerical factors); 
in Sec. IVI1 we will consider verifications of macroscopic 
quantum entanglement between test masses in GW de- 
tectors as a test of gravity decoherece; in Sec. lVIIl we will 
summarize our main results. In the Appendix, we will 
present mathematical details for solving integral equa- 
tions that we encounter in obtaining the optimal verifi- 
cation scheme. 



II. MODEL AND EQUATIONS OF MOTION 

In this section, we will present a mathematical descrip- 
tion of the system model, as shown schematically in the 
upper left panel of Fig. [SJ The oscillator position is lin- 
early coupled to coherent optical fields through radiation 
pressure. Meanwhile, information of the oscillator posi- 
tion flows into the outgoing optical fields continuously. 
This models a measurement process in an optomechan- 
ical system without a cavity or with a large-bandwidth 



cavity. The corresponding Heisenberg equations, valid 
for both preparation and verification stages, are formally 
identical to classical equations of motion except for that 
all quantities are Heisenberg operators. The oscillator 
position x and momentum p satisfy the following equa- 
tions: 



x(t) = p(t)/m, 



(12) 



P (t) = -2 7m p(t) - mu) 2 m x{t) +aa 1 (t)+£ F (t). (13) 

Here a d\ corresponds to the quantum-radiation-pressure 
noise or so-called back-action noise; a = {km SI 2 ,) 1 / 2 = 
(8 Iq ujq h/c 2 ) 1 / 2 is the coupling constant between the 
oscillator and optical fields with Iq denoting the op- 
tical power and Vl q quantifying the characteristic fre- 
quency of measurement strength. We have included the 
fluctuation-dissipation mechanism of the mechanical os- 
cillator by introducing the mechanical damping rate "f m 
and classical- force noise i-C, the Brownian thermal 
noise. In the Markovian limit, the correlation function 
for £f is given by 1 



(ZF(t)ZF(t')) s ym = Sf5(t-t')/2 



(14) 



where Sfi 



4m"f m kBTo = 2hmfl 2 7 and we have defined 



a characteristic frequency £If for the thermal noise. 

The amplitude and phase quadratures of ingoing op- 
tical fields ai,2 and of outgoing optical fields 61,2 satisfy 
the following input-output relations: 



h{t) = y/rjni(t) + x/l-r?ai(i), 



b 2 (t) = 0?n 2 (t) + v 7 !-?/ |S a (t) + jr£(t) + 



(15) 

1 . 
(16) 



Here h\ t 2 originate from nonunity quantum efficiency of 
the photodetcctor for r\ > 0. In the paraxial and narrow- 
band approximation, ai.2 are related to the electrical- 
field strength at the central frequency ujq by [52l454| : 



Elt) = — {la + aim coswo* + a 2(i) smcj i} 

\ Sc J 

(17) 

with a denoting the classical amplitude and S standing 
for the effective cross-section area of the laser beam. A 
similar relation also holds for the outgoing fields 61,2- 
In addition, they satisfy [ai(i), a^t 1 )] = [fei(i), 62^')] = 
iS(t — <'). Their correlation functions read 



(«,(*)%■((')), 



6 l3 e ±2 i6(t-t')/2, (i,j = l,2) (18) 



Here ( } S ym stands for a symmetrized ensemble average. For a 
system characterized by a density matrix p, it is defined as 

(61 (t) 02(t'))sym = Tr {[oi(t)6 2 (t') + d 2 {t')o 1 (t)]p} /2 . 
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FIG. 5: A schematic plot of the system (upper left panel) and the corresponding spacetime diagram (right panel) showing 
the timeline of the proposed MQM experiment (see Sec. IIII Al for detailed explanations). In this schematic plot, the oscillator 
position is denoted by x which is coupled to the optical fields through radiation pressure. The ingoing and outgoing optical 
fields are denoted by ai,2 and 61,2 with subscripts I, 2 for the amplitude and phase quadratures, respectively. In the spacetime 
diagram, the world line of the oscillator is shown by the middle vertical line. For clarity, ingoing and outgoing optical fields are 
represented by the left and right regions on the different sides of the oscillator world line, even though in reality, optical fields 
escape from the same side as where they enter. We show light rays during preparation and verification stages in red and blue. 
In between, the yellow shaded region describes the evolution stage with light turned off for a duration of te- The conditional 
variance of the oscillator motion is represented by the shaded region alongside the central vertical line (not drawn to the same 
scale as light propagation) . At the beginning of preparation, the conditional variance is dominated by that of the initial state 
(orange). After a transient, it is determined by incoming radiations and measurements. Right after state preparation, we show 
the expected growth of the conditional variance due to thermal noise alone, and ignoring the effect of back-action noise which 
is evaded during the verification process. The verification stage lasts for a duration of tv, and it is shorter than tf after which 
the oscillator will be dominated by thermalization. 



where q denotes the squeezing factor (q = for a vacuum- 
state input) with "+" for the amplitude quadrature and 
"— " for the phase quadrature. Correspondingly, the cor- 
relation function for the back-action noise a a\ is simply 

(aa 1 (t)aa 1 (t')) sym = S3 A S(t-^)/2 (19) 

with SfS A = e 2q hmil 2 q . In Eq. ([16]). £ x is the sensing 
noise. One example is the internal thermal noise, and it 
is defined as the difference between the center of mass 
motion and the surface motion of the oscillator which is 
actually being measured. In the Markovian approxima- 
tion, it has the following correlation function: 

(i(t)i(t'))sym = Sf^-0/2 (20) 

where <S* h = h/(mfl x ) and we introduce a characteristic 
frequency Q x for the sensing noise. 

Note that the fi g , Hp, an d Q x that we have introduced 
are also the frequencies at which the back-action noise, 
thermal noise and sensing noise intersect the SQL [cf. 
Eq. (JTJ], respectively. They are identical to what were 
introduced in Ref . [35| . For conveniences of later discus- 
sions, we introduce the following dimensionless ratios: 

Cf = n F /n q , ( x = n q /n x . (21) 

In addition, we define two characteristic timescalcs for 
the measurement and thermal-noise strength as 

r q = I/O, , t f = i/n F . (22) 



III. OUTLINE OF THE EXPERIMENT WITH 
ORDER-OF-MAGNITUDE ESTIMATE 

In this section, we will describe in details the timeline 
of a plausible MQM experiment (subsection IIII A[) and 
provide order-of-magnitude estimates of the conditional 
variance of the prepared quantum state, the evolution of 
the prepared quantum state, and the verification accu- 
racy in the free- mass regime (subsections IIII B[ IIII CI and 
IIII Dp . This will provide qualitatively the requirements 
on the noise level for the success of a MQM experiment. 
We will give more rigorous treatments in the next section. 



A. Timeline of proposed experiment 

We have sketched a space-time diagram for the pro- 
posed MQM experiment in the right panel of Fig. [5] — 
with time going upward, therefore we start from the bot- 
tom of the figure. 

Lock Acquisition. At the beginning, the mechanical os- 
cillator is in a highly mixed state, so are the optical 
fields. Therefore, the first step is to "acquire lock" of the 
measurement device, and reach a steady-state operation 
mode, during which several r q will have elapsed. From 
this time and on, initial-state information will have been 
forgotten (propagating outward within the green strip), 
and the state of the oscillator will be determined by the 
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driving fields, including the classical-force noise and sens- 
ing noise, as well as the quantum noise. This will be the 
start of the state-preparation stage (region above the 45° 
green strip). 

State Preparation. This stage is a steady-state opera- 
tion of the measurement device. The quantum state of 
the oscillator is collapsed continuously due to homodyne 
readouts of the photocurrent. At any instant during state 
preparation, based on the measured history of the pho- 
tocurrent (mostly on data within several times r q to the 
past oft), the conditional expectation (a; cond , p cond ) for 
the oscillator position x and momentum p can be con- 
structed. The second moments, dcscribable by the co- 
variance matrix between position and momentum, which 
consists of V; c ° nd , V"™ nd and V™ nd , can be calculated 
from the noise model of the measurement device — they, 
together with x cond and p cond 7 folly determine the quan- 
tum state, i.e., the Wigner function of the oscillator at 
any instant [cf. Eq. @]. For a Gaussian steady state, the 
construction of (x cond , p cond ) and conditional covariancc 
matrix from the history of the photocurrent can be ac- 
complished most easily using Wiener Filtering, as shown 
in Ref. (H- 

The preparation stage terminates at t = 0, when 
^ x cond^ pcond-j an( j covar j ance matrix will be determined 

by data from several — r q up to as shown by the red 
strip. 

State Evolution. If we want to investigate the quantum 
dynamics of the oscillator and study various decoherence 
effects, we can delay the verification process and allow the 
oscillator to freely evolve with the interaction light turned 
off (represented by the yellow strip). During this period, 
the thermal noise will induce diffusions of the oscillator 
position and momentum, thus increasing the conditional 
variance as shown schematically by broadening of the 
shaded region alongside the oscillator world line. If there 
were any additional decoherence effect, the variance will 
grow faster than the case with the thermal decoherence 
alone. A follow-up verification allows us to check whether 
additional decoherence mechanisms, such as the gr avity 
decoherence conjectured by Diosi [49| and Penrose [501 ] . 
exist or not. 

State Verification. After the evolution stage, the verifi- 
cation stage starts (represented by blue strip). We in- 
tentionally use different colors to label the preparation 
light and verification light — symbolizing the fact that 
in principle, a different observer (verifier) could perform 
the verification process, and verify the quantum state 
by him/herself. The only knowledge from the preparer 
would be the conditional expectation x cond and p cond 
if all noise sources are Markovian. The verifier uses 
a time-dependent homodyne detection and collects the 
data from measuring the photocurrents. The verification 
process lasts for a timescale of ry between the charac- 
teristic measurement timescale r q and the thermal deco- 
herence timescale Tp, after which diffusions of x and p 
in the phase space become much larger than the Heisen- 
berg limit. Based upon the measurement data, the veri- 



fier can construct an integral estimator for one particular 
mechanical quadrature [cf. Eq. [fj]. 

The above three stages have to be repeated for many 
times before enough data are collected to build up statis- 
tics. After finishing the experiment, the verifier will ob- 
tain a reconstructed quantum state of the mechanical 
oscillator, and then can proceed to compare with the pre- 
parer and interpret the results. 



B. Order-of- magnitude estimate of the conditional 
variance 

In this and the following two subsections, we will 
provide ordcr-of-magnitudc estimates for a three-staged 
MQM experiment including preparation, evolution and 
verification stages. This gives us physical insights into 
different timescales involved in a MQM experiment and 
also the qualitative requirements for an experimental re- 
alization. We will justify those estimates based upon 
more careful treatments in the next several sections. 

Based upon the measurement data from several — r q to 
0, one can construct a conditional quantum state for the 
mechanical oscillator. Suppose that the phase quadra- 
ture of the outgoing fields is being measured and the 
photodetection is ideal with r} = 0. Given a measure- 
ment timescale of r (measuring from — r to 0), variances 
for the oscillator position and momentum at t = in the 
free-mass regime with tu m — > are approximately equal 
to [cf. Eqs. |Q2), 03), US) and (USD] 

fe 2 (0)-5* ot /r + t 3 S*£7to 2 ~ NiNj8x 2 q , (23) 
5p 2 (0)~m 2 Sl ot /T 3 + tS^ ~ NiNjSp 2 . (24) 

Here S^ ot = Sf A + Sf [cf. Eqs. and (QII ] and 
5* ot = S x h + with S s x h denoting the shot noise due to 
a 2 [cf. Eqs. (fTo) and ([2"U|) ]; we have defined 

N x = l + 2( 2 , Np = 1 + 2(1, (25) 

while 

5x 2 q = h/(2mQ q ) , Sp 2 q = hmn q /2. (26) 

The optimal measurement timescale is given by r ~ r q . 
Purity of the prepared conditional quantum state at t = 
is approximately equal to [cf. Eq. (|3J)] 

17(0) ~ hx(0) Sp(0) ~ N x Np. (27) 

If classical noises are low, namely, N x ~ Np ~ 1 , the con- 
ditional quantum state will be pure with U(0) ~ 1. For 
future GW detectors such as AdvLIGO, both ( x and (p 
will be around 0.1, and such a low classical- noise budget 
clearly allows us to prepare nearly pure quantum states 
of the macroscopic test masses. 
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p/(mu; m ) 

t = o- 




Thc growth in the uncertainty ellipse will simply be 

U th (r E ) ~ ^fe th (r B ),5p th (T B ) ~ C 2 F (n q T E ) 2 = {t e /t f )\ 

(32) 

When t e > Tf, U th (r E ) > 1 and the conditional quan- 
tum state will be dominated by thcrmalization. 

If there were any additional dccoherence effect, the 
growth in the uncertainty will be much larger than what 
has been estimated here. A subsequent verification stage 
can serve as a check. 



FIG. 6: (Color online) Rotation and diffusion of a highly 
position-squeezed conditional quantum state prepared by a 
strong measurement with Q q 3> c^m- The initial- momentum 
uncertainty will contribute an uncertainty in the position 
comparable to the initial position uncertainty when the evo- 
lution duration te ~ r q - 



C. Order-of-magnitude estimate of state evolution 

During the evolution stage, the uncertainty ellipse of 
the conditional quantum state will rotate at the mechani- 
cal frequency in the phase space, and meanwhile there is a 
growth in the uncertainty due to thermal decoherence as 
shown schematically in Fig. [5J Given a strong measure- 
ment, the variance of the resulting conditional quantum 
state in position fe 2 (0) will be approximately equal to 
8x 2 as shown explicitly in Eq. ([2"B")) with N x ,Np ~ 1. 
It is much smaller than the zero-point uncertainty of 
an u m oscillator, which is given by h/(2moj m ). There- 
fore, the conditional quantum state of the oscillator is 
highly squeezed in position. The position uncertainty 
contributed by the initial-momentum will be comparable 
to that of the initial-position uncertainty after a evolu- 
tion duration of T q . This can be directly seen from an 
order-of-magnitude estimate. In the free-mass regime, 



x(t) ~ .t(0) + ^-t. 

m 



(28) 



For an evolution duration of t e , the corresponding vari- 
ance in position is 



6x 2 (t e ) ~ fa 2 (0) + 



5p 2 (0) 



8x 2 {0)[l + (n q T E ) 2 



(29) 

The contribution from the initial-momentum uncertainty 
(the second term) will become important when fl q r E ~ 1 
or equivalently t e ~ r q . 

Apart from a rotation, the uncertainty ellipse will also 
grow due to thermal decoherence. Variances in the posi- 
tion and momentum contributed by thermal decoherence 
are approximately given by [cf. Eqs. (|T2l and ([T3"|) ] 



T E Sf/m 2 = C F (n q T B rSx' q , (30) 
t e Sf = ( F (n q T E )6p 2 q . (31) 



D. Order-of-magnitude estimate of the verification 
accuracy 

To verify the prepared conditional quantum state, the 
oscillator position needs to be measured for a finite du- 
ration to obtain information about x(0) and p(0) [cf. Eq. 
([5]) and ©] or about x(t e ) and p{t e ) if the evolution 
stage is inserted. In order for an entire state character- 
ization to be possible, one might then expect that an 
oscillation period must pass, and during this period, the 
thermal noise should cause an insignificant diffusion of 
the oscillator momentum c omp ared with its zero-point 
uncertainty, which requires [40j 



< 



(33) 



with Q m = cu m /(2"fm) denoting the mechanical qual- 
ity factor. This requirement is unnecessary if the ini- 
tial quantum state is prepared by a strong measurement. 
As we have mentioned in the previous subsection, the 
resulting condition quantum state is highly squeezed in 
position and the initial-momentum uncertainty will make 
a significant contribution to the uncertainty in position 
after r > r q . This means, depending on the particular 
strategy, one can extract x and p below the levels of Sx q 
and 8p q , respectively, as long as one is able to measure os- 
cillator position with an accuracy better than 8x q , within 
a timcscale of several r q . This is certainly possible if the 
measurement-induced back action is evaded. 

To evade the measurement-induced back action, one 
notices the fact that the amplitude quadrature b\ con- 
tains d\ which is responsible for the back action, and 
meanwhile the phase quadrature 62 contains the infor- 
mation of oscillator position, part of which is contributed 
by the back action [cf. Eq. ([TS]) - (|T5|) ] . Therefore, if we 
measure particular combinations of b\ and hi at different 
times, by summing up those measurements, we will be 
able to cancel the back action and obtain a back-action- 
evading (BAE) estimator for a given mechanical quadra- 
ture. Such cancelation mechanism is only limited by the 
readout loss (77 ^ 0), which introduces uncorrelated vac- 
uum fluctuations. 

We can make an order-of-magnitude estimate to show 
that a sub-Hcisenberg accuracy can be indeed achieved. 
With the BAE technique, the force noise that limits the 
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verification accuracy will only contain the thermal-noise 
part. Similar to Eqs. {23} and (JUJ) but with replaced 
by Sfi, the variances in position and momentum during 
the verification stage are simply 



Sxy 
Sp 2 v 



m 2 S t ° t /T 3 + T Sf ~ 



(34) 
(35) 



Here the optimal verification timescale would be ry ~ 

— 1/2 

C F T q and T q < Ty < Tp. A summarizing figure of 
merit for the verification accuracy is approximately given 
by 



j radd I 

U BAE 



-SxySpv 
h 



(36) 



A sub-Heisenberg accuracy can be achieved when < 1. 
Note that this error can be arbitrarily small by lowering 
£f indefinitely, i.e., a very strong measurement. If phase- 
squeezed light is injected during the verification stage, we 
would have 



T j add I 

U BAE 



2U/2, 




Increasing the squeezing factor always improves our ver- 
ification sensitivity, with a limit of 



^lini IbAE ~ FIf/Qx — Cx (f, 



(38) 



which can be much lower than unity in the case of future 
GW detectors or any low-noise measurement device. 

Had we not evaded the back-action noise, we would 
have y/Np in the place of which means SxySpv would 
be Heisenberg-limited — unless different squeezing fac- 
tors are assumed. For low squeezing (i.e., e ±2q larger 
than both C, x and Cf)> we need phase squeezing for x 
observation, amplitude squeezing for p observation, with 



the most sensitive frequency band (f2 ~ fl q ) of the mea- 
surement device, not only can we prepare a nearly pure 
quantum state, but also can we make a sub-Heisenberg 
tomography of the prepared state. In this and following 
sections, we will provide more rigorous treatments di- 
rectly by analyzing the detailed dynamics of the system. 



A. The conditional quantum state obtained from 
Wiener filtering 

The rigorous mathematical treatment of state prepa- 
ration has been given in Ref. [35j . The main idea is 
to treat the conditional quantum state preparation as a 
classical filtering problem, which is justified by the fact 
that the outgoing optical quadratures at different 
times commute with each other, the same as a classi- 
cal random process. For such a Gaussian linear system, 
the Wiener filter, satisfying the minimum mean-square 
error criterion, allows us to obtain an optimal estimate 
for the quantum state of the oscillator, i.e., the condi- 
tional quantum state. Based upon the measurement data 
y(t) (t < 0), conditional means for oscillator position and 
momentum at t = can be constructed as [cf. Eq. (14) 
of Ref. M] 



x cond (0) = (x(0)Y 



p cond (0) = (p(0)) 



cond 



dtK x (-t)y(t), (41) 



dtK p (-t)y(t). (42) 



Here K x and K p are causal Wiener filters. The covariance 
matrix is given by [cf. Eq. (15) of Ref. HH ] 

v^(o) = (0,(0)0,(0))™;;? - ( Ol (o)r nd (o,(o)r nd , 

(43) 

where i,j = 1,2 and 01,62 denote respectively. In 
the free-mass regime, we showed that [cf. Eq. (52)-(54) 
in Ref. M]: 



ithout BAE 



(39) 



which is a significant factor (1/Cf) worse than the BAE 
scheme. Even though there exists an optimal squeezing 
factor that this scheme can apply and yields 



7-7-add 1 

U opt I without BAE 



(40) 



yet it is still worse than the limiting situation of the BAE 
scheme [cf. Eq. (|3"8"])] by a factor of l/Q F (> 1). 



IV. 



THE CONDITIONAL QUANTUM STATE 
AND ITS EVOLUTION 



The previous order-of-magnitude estimates provide us 
a qualitative picture of a MQM experiment, especially 
in the free-mass regime where future GW detectors are 
operating. As long as and Q x are smaller than unity, 
namely, the classical noise goes below the SQL around 



V 



cond 



(0) 



N~pN, 
Nl 



~V28x 2 q 



NpN x 2 h/2 



(44) 



, F mh/2 N F NiV28p 2 

With conditional means and variances, the Wigner func- 
tion or equivalently the conditional quantum state is 
uniquely defined [cf. Eq. @]. Correspondingly, purity 
of the conditional quantum state is quantified by 



U(0) = Wdet V cond (0) = N X N F . 



(45) 



This simply justifies the order-of-magnitude result pre- 
sented in Eq. (|27|) . 



B. Evolution of the conditional quantum state 

In the following discussions, we will analyze how such 
a conditional quantum state evolves during the evolu- 
tion stage. On the one hand, this confirms the quali- 
tative results presented in the subsection MI CI On the 
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other hand, it provides a quantitative understanding of 
the timescale for the later verification stage. 

The equations of motion for the oscillator during the 
evolution stage are given by Eqs. (fT2|) and (fT3)) except 
that there is no radiation pressure, for the light is turned 
off 2 . For simplicity and also a consideration of the case 
in a realistic experiment, we will assume an oscillator 
with a high quality factor, i.e., w m 3> 7 m . Within a 
timescale much shorter than l/j m , the oscillator can be 
well-approximated as a free oscillator. Correspondingly, 
the analytical solution to oscillator position reads 



x(t)=x q (t) + dt'G x (t-t')Z F (t'). (46) 
Jo 



Here the free quantum oscillation x q (t) of the oscillator is 
given by Eq. ([5]). We have defined the Green's function 



G x (t) = e(t) 



sin(u m t) 
muj m 



(47) 



with O(i) denoting the Heaviside function. 

Given an evolution duration of t e , from Eq. (|14[) and 
Eq. (|46|) the corresponding covariance matrix evolves as 



Rl V cond 

Cth 



8 mru) 



2,. ,3 



(0) R 

2 $ - 

2 mw. 



sin 2 $ 
,sin 2 3> 



sin 2 $ 



2 moj. 

2 (2$ + sin 2$) 
(48) 



where $ 
by 



w m T£ and the rotation matrix R$ is given 



series of Up to the leading order in we obtain 



V xx {t e ) = V^ + 



4<5.t 2 Sx^ 

q Tycondo ~- i ^ T/cond/o 



2Sx 2 X 



V xp (t e ) = V™ d + 



3 



T/condo 



v pp (t e ) = v p c ° nd + 2Sp 2 q ( F n q r E 



(50) 



(51) 
(52) 



denoting the elements of V cond (0). Up to 
the leading order in 0, q r E , the uncertainty product of the 



with V con ' 



resulting quantum state is 



U(t e ) = -VdetV(r E ) 



[7(0) 



T/cond 

~tf {te/tf)2 (53) 



with t e defined in Eq. (|22|) . The second term is con- 
tributed by the thermal decoherence and can be viewed 
as U th (r E ). Those formulas recover the results in Eqs. 
(f2T)f - ([3"!?)) but with precise numerical factors. As we can 
conclude from Eq. (|53[) . in order for a sub-Hcisenberg 
tomography to be possible, the later verification stage 
should finish within a timescale of tf after which the 
contribution from the thermal noise gives II th (t e ) ~ 1. 



V. STATE VERIFICATION IN THE PRESENCE 
OF MARKOVIAN NOISES 

In this section, we will treat the followup state verifi- 
cation stage with Markovian noises in details. This can 
justify the order-of-magnitudc estimate we have done in 
the subsection IIII Dl In addition, we will show explicitly 
how to construct the optimal verification scheme that 
gives a sub-Heisenberg accuracy. 



R 



* — 



(m w 



cos$ — 7noj m sin$ 
~ 1 sin$ cos<3> 



(49) 



A time-dependent homodyne detection and 
back-action-evading (BAE) 



The first term in Eq. (|4"5)l represents a rotation of the 
covariance matrix V cond (0) due to the free quantum os- 
cillation of the oscillator; the second term is contributed 
by thermal decoherence which causes an increase in the 
uncertainty. 

In the free-mass regime and the case of uj m T E <C 1, 
elements of the covariance matrix can be expanded as 



2 Were the light turned on, the back action can still be evaded as 
long as one measures the amplitude quadrature ai during this 
period and take them into account during data processing. Since 
no information of the oscillator position (contained in the phase 
quadrature of outgoing light) is collected, this is equivalent to 
the case with light turned off. 



In this subsection, we will analyze the time-dependent 
homodyne detection which enables us to probe mechan- 
ical quadratures. We will further show how the BAE 
scheme can be constructed. The BAE scheme is optimal 
only when there is no readout loss (rj = 0). We will con- 
sider more general situations and derive the correspond- 
ing optimal verification scheme in the next subsection. 

The equations of motion for the oscillator during the 
verification stage (t > t e ) are given by Eqs. (|12[) and 
([T3l) . The corresponding solution to oscillator position 
is different from Eq. (|46[) due to the presence of the 
back-action noise which starts to act on the oscillator at 
t = t e . Specifically, it reads 

P oo 

x(t)=x q (t)+ / dt / G x (t-t / )[«ai(* / )+Cp(* / )]- (54) 
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FIG. 7: (Color online) A schematic plot of time-dependent 
homodyne detection. The phase modulation of the local os- 
cillator light varies in time. 



Here the free quantum oscillation x q (t) is the signal that 
we seek to probe during the verification stage. For op- 
tical quadratures, the equations of motion are given by 
Eqs. (fT5)) and (fT5|) . From those equations, we notice that 
among the outgoing fields: b\ is pure noise, while b 2 con- 
tains both signal x q (t) and noise. In order to highlight 
this, we rewrite &1.2 as 



61 (t) = y/1 - rjntit) + y/rj&i(t) = Shit) , (55) 

Sa(t) = Sk{t) + y/T^i{a/h)x 9 (t) (56) 
with [cf. Eq. ([53])] 

6b 2 (t) = y/rjM*) + V 7 ! - v{a2(t) + f 

At' G x (t-t') [a oi(tO +&■(«')]}• 

(57) 



In this way, we can directly see that d\ which causes the 
back action is contained in both the amplitude quadra- 
ture 61 and the phase quadrature 62- Therefore, by 
measuring an appropriate combination of the two out- 
put quadratures, we will be able to remove effects of the 
back-action noise that is imposed onto the oscillator dur- 
ing the verification process at t > te- Searching for such 
an optimal combination is the main issue to be addressed 
in this section. 

As mentioned in the introduction part, to probe me- 
chanical quadratures and their distributions, a time- 
dependent homodyne detection needs to be applied [cf. 
Eq. (JBJ)]. Specifically, the outgoing optical field 



B ou t(t) = b\(t) cos u} t + b 2 (t) sinw i 



(58) 



at t > te is mixed with a strong local-oscillator light 
L(t) whose phase angle <p os is time-dependent as shown 
schematically in Fig. namely, 



L(t) = L cos[w 1 - 4>os(t)} 



(59) 



with Lq a time- independent constant. Through a low- 
pass filtering (with a bandwidth much smaller than ujq) 
of the beating signal, the resulting photocurrent is 



i(t) oc 2B out (t)L{t) 

= L bi (t) cos (j) os (t) + L 62 (i) sin (f> os (t) , (60) 



where the overline means averaging over many optical- 
oscillation periods. Note that Heisenberg operators for 
the photocurrent at different time commute with each 
other, i.e., 



[i(t)Mt')] = o. 



(61) 



and are therefore simultaneously measurable, as obviously 
expected. Based on the measurement results of i(t) from 
te to Ti n t, we can construct the following weighted quan- 
tity Y with a weight function W(t): 

/•Tint 

Y = / e(t-T E )W{t)i{t)dt = {gi\b 1 ) + (g 2 \b 2 ) . (62) 
Jo 

Here the Heaviside function Q(t — te) manifests the fact 
that the verification stage starts at t = te and we have 
introduced the scalar product of two vectors \A) and \B) 
in the £ 2 [0,Ti n t] space as the following: 



(A\B) 



A(t)B(t)dt. 



(63) 



Besides, we have defined filtering functions g\ and g 2 as 

gi (t) = e(t-T B )W(t)cos4> os (t), (64) 
g 2 (t) = Q(t - T E )W(t) sin0 os (i). (65) 

Since all the data can in principle be digitalized and 
stored in hardwares, the weight function Wit) can be 
realized digitally during data processing. In addition, an 
overall re-scaling of gi, 2 it) — > Co gi, 2 {t) with Co a time- 
independent constant docs not affect the verification per- 
formance, and that there are multiple ways of achieving 
a particular set of gi^it), by adjusting the phase </> os (i) 
of the local oscillator and the weight function W(t). 

In light of Eqs. (|55|) - (|57| . we decompose the weighted 
quantity Y [cf. Eq. (|62"]) ] as a signal Y s and a noise part 
SY, namely, 

Y = Y s + SY. (66) 

They are given by 

Y s = y/l-r)ia/h) ig 2 \x g ), 

SY = (g^Sh) + iffilSk)- (67) 

Since an overall normalization of gi, 2 will not affect the 
signal-to-noise ratio as mentioned, we can impose, math- 
ematically, that 



(52I/1) =cosC, (521/2) = sinC 



(68) 



with 



flit) = cosw m t, / 2 (t) = iCL q /u m )sinu} m t (69) 

in the coordinate representation. The signal part can 
then be rewritten as 



Y s = y/l - i) ia /K)5x q [x cos ( + po sin £] , (70) 
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where we have introduced normalized oscillator position 
and momentum as Xq = x(te)/Sx q and po = p(rE)/Sp q . 
In such a way, a mechanical quadrature of Xq will be 
probed [cf. Eq. (J2)]. For the noise part, more explicitly, 
we have [cf. Eqs. ([55]) - ([57]) ] 

$Y = {gAVvni + V 1 - T l a i) + {g2\Vv n 2 + V 1 ^^^) 

+ y/T^ (a/h)[(g 2 \G x \i F ) + (g 2 \i x )} , (71) 

where the integration with G x (t—t') has been augmented 
into applying a linear operator G x in the £ 2 [0, Tj nt ] 
space. In the above equation, terms on the first line 
is the shot noise, and the term on the second line is the 
back-action noise, while terms on the third line are the 
classical-force and sensing noises. 

The optimal gi (t) and g 2 (t) that give a sub-Heisenbcrg 
accuracy for each quadrature will be rigorously derived 
for general situations in the next section. If ai and a 2 
arc uncorrelatcd and there is no readout loss with rj = 0, 
an optimal choice for g\ would be obvious to cancel the 
entire contribution from the back-action noise term (pro- 
portional to fii). This is equivalent to impose, mathemat- 
ically, that 

(g 1 \a 1 ) + (a 2 /h)(g 2 \G x \a 1 ) = (72) 

or 

\g 1 ) + (a 2 /h)Gf>\g 2 ) = 0, (73) 

where G| J is the adjoint of G^. Physically, this cor- 
responds to bringing in a piece of shot noise (gi\ai) to 
cancel the back-action noise (a 2 /h){g 2 \G x \a\) — there- 
fore achieving a only shot-noisc-limitcd measurement. In 
the coordinate representation, Eq. (|73p can be written 
out more explicitly as 

/Tint 
dt'G x (t' -t)g 2 (t') = 0, (74) 

which agrees exactly with the variational-type BAE mea- 
surement scheme first investigated by Vyatchanin et al. 
[42| . It is suitable for detecting signals with known arrival 
time. For stationary signals, one would prefer frequency- 
domain variational techniques proposed by Kimble et 
al. [13] , which evades the back-action noise for all possible 
signals as long as they are Gaussian and stationary. 

As realized by Kimble et al. [52} in their frequency- 
domain treatment, when the readout loss is significant 
(large rf) and when the back-action noise is strong (large 
a), the variational approach becomes less effective, be- 
cause in such a case, the magnitude of g\ required to 
bring enough a\ to cancel the back-action noise would 
also introduce significant noise hi [cf. Eq. ([71])]. This 
reasoning apparently leads to a trade-off between the 
need to evade back action and the need to minimize loss- 
induced shot noise — such an optimization will be made 
in the next section. 



B. Optimal verification scheme and covariance 
matrix for the added noise: formal derivation 

Imposing the BAE condition [cf. Eq. ([741 ] does not 
specify the shape of g 2 , nor does Eq. ([55]) . and we have 
further freedom in choosing g 2 that minimizes the noise 
in measuring a particular quadrature of Xq . In addition, 
in the presence of readout loss with r\ ^ 0, totally evad- 
ing back action is not the obvious optimum as mentioned. 
Therefore we need to optimize g\ and g 2 simultaneously. 
In this section, we first carry out this procedure formally, 
and apply to the Markovian-noise budget in the next sub- 
section. 

The total ig-referred noise in the weighted estimator 
Y can be written as [cf. Eqs. ([70]) and ([7T]) ] 

^= ( i4W ^ )syro 

= ^^ q t^\CM, (75) 

where correlation functions Cy- among the noises are the 
following: 

C ij (t,t') = (^(t)5S i (t , ))sy I n ! (i,j = l,2). (76) 

The optimal gi, 2 (t) that minimize a 2 can be obtained 
through the standard constraint variational method. For 
this, we define an effective functional as 

J ctf = (l-77)(^/4)tr 2 [ 5 i )2 ]- i ui(/i|32)-M2(/2|52) 

= lj2(9i\ C ij\9j) ~ Gui/i+M2/2|<?2), (77) 

id 

where fj,i and fj, 2 are Lagrange multipliers due to the 
normalization constraints in Eq. (|68[) . Requiring the 
functional derivative of J e R with respect to g\ and g 2 
equal to zero, we obtain 

Cii|3i) + Ci a |5a) = 0, (78) 
C 2 i|3i) + C 22 | ff2 ) = |/ii/i + M 2 / 2 ). (79) 

Here should be viewed as operators in the £ 2 [0, Tint] 
space. This leads to formal solutions to g\ 2 , namely, 

= -CtfCulsa) , (80) 
\g 2 ) = M| Ml A+/i 2 / 2 ), (81) 

where we have defined 

M= [C^-C^C^C!,]" 1 . (82) 

Re-imposing Eqs. (|68[) . those unknown Lagrange multi- 
pliers [i\ %2 can be solved, which are related to C by 



(A|M|A) 


(/i|M|/ 2 )" 








cos£ 


(/ 2 |M|/i) 


(/ 2 |M|/ 2 ) 




^ 2 _ 




sin £ 



(83) 
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Correspondingly, the minimum cr^ in has the following 
quadratic form: 



^min = I C0S C sinC]V 



add 
norm 



cos( 
sin £ 



(84) 



Here normalized V*o d m 



is a 2 x 2 covariance matrix, and 



it is given by 



V 



add 



(i - v )n q 



(/i|M|/i) (/i|M|/ 2 ) 
(/a|M|/i) (/ 2 |M|/ 2 ) 



(85) 



It relates to the initial definition of the covariance matrix 
for the added verification noise [cf. Eq. {5)] simply by 



V add = Diag[fe 3 , fcjV^Diagp*,, S Pq 



(86) 



Due to the linearity in Eqs. ((79)) and {53"). the opti- 
mal <7i.2 for a given quadrature ( can also be rewritten 
formally as 



51,2 = 9l,2 C0S C + 9l,2 sin C : 



(87) 



with p 12 = gl i2 (0) and g{ 2 = g\ 2 (iz/2). Such 
(■-dependence of gi t2 manifests the fact that a sub- 
Heisenberg tomography requires different filtering func- 
tions, or equivalcntly different measurement setups, for 
different quadratures. 



C. Optimal verification scheme with Markovian 
noise 

Given Makovian noises, the corresponding correlation 
functions for the output noise dbi can be written out ex- 
plicitly as [cf. Eqs. (H), {H), {2U), and {76)] 



Cn(t,fl ^ + V )e2 V -0, 



C 12 (t,t') = C 21 (t',t) = (l-^—^G^t'-t), (89) 



A 2 



a 4 /e 2c ' 



c 22 (t,t') = ^s(t-t') + (i- v y^ r—+c 2 F 



dt x G x {t-t x )G x {t' -t x ), 



(90) 



with A = ^2[rj + (1 - 7y)(e- 2 9 + 2C 2 )]. Plugging these 
Cjj into Eq. ([8"U|) and (|8T|) . we can obtain the equations 
for the optimal filtering functions g x and g 2 . Specifically, 
for g x , we have [cf. Eq. (|8"0) ] 

fl - nk 29 a 2 /" Tint 

(91) 

For ^2, by writing out M explicitly, it gives [cf. Eq. {81)] 



A 2 „„a 4 



4—' ' " fi 2 JJ0 

= Mi/i(i)+/"2/ 2 (t), 



dt'dt x G x {t - t x )G x (t' - h)g 2 (t') 

(92) 





1.0 




0.5 


cf 






0.0 












0.5 




1.0 











I \ 
/ \ 




-"*' P 













ttqt 



10 



15 



FIG. 8: (Color online) Optimal filtering functions gi (solid 
curve) and gi (dashed curve) in the presence of Markovian 
noises. We have assumed Q. q /2n = 100 Hz, C, x = C,f = 0.2, 
r\ — 0.01 and vacuum input (q — 0). For clarity, the origin of 
the time axis has been shifted from te to 0. 



where we have introduce £' F , which is given by 



77(1 — rf)e 



2<l 



2[r/+ (1 -r))e 2 i] 



+ (1 - v)C 2 f 



1/2 



2 +Cf 



1/2 



(93) 

and it is equal to (f for no readout loss. Although here g\ 
is still defined from g 2 , the optimal verification strategy 
does not totally evade the back action, as is manifested in 
the term proportional to r] inside the bracket of Eq. (|93p . 
In the limit of no readout loss with rj = 0, it is identical to 
the BAE condition in Eq. {74) . Typically, we have 1% 



readout loss X) = 0.01, squeezing e 2q = 10 and Cf = 0.2, 
this readout loss will only shift O by 6%, which is neg- 
ligible. However, if the thermal noise further decreases 
and/or the measurement strength increases, the effect of 
readout loss will become significant, entering in a similar 
way as the frequency-domain variational measurement 
proposed by Kimble et al. [HJ . 

The above integral equations for optimal g\ and g 2 
can be solved analytically as elaborated in the Appendix 
[Cl which in turn gives M and the corresponding V add 
[cf. Eqs. (|52") and ([55]) ]. In the free-mass regime with 
fl q uj m , closed forms for optimal g\ and g 2 can be 
obtained, which, in terms of g 1 2 [cf. Eq. (|87|) 
given by 



are 



</f = .9i| c =o = (n q /x)e- n " xt sinn q xt; (94) 
fff = </i|c= f = -V2Q q e-n^sin (fi, X t + J) , (95) 

and 



02 =.92lc=o = 20 g xe "^cosQ^t 
.9 2 P =fl2k=f = 2V2fi ?X 2 e 



(96) 



* = 2V2 r! oX 2 e-°' xt sin ( n qX t - - ) , (97) 



with x = [C^/A] 1 / 2 . The corresponding verification 
timcscale is set by ry = (xfi g ) _1 and T g < Ty < t/t. To 
illustrate the behavior of the optimal filtering functions, 
we show g\' 2 P in Fig. [5] As we can see, the verification 
process finishes after several r q , i.e., in a timescale of ry. 
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-1.5-1.0-0.5 0.0 0.5 1.0 1.5 

x/5x q 

FIG. 9: (Color online) The uncertainty ellipse for the added 
verification noise in the presence of Markovian noises. We 
assume ( x — £p = 0.2, vacuum input (Dashed curve), £ a = 
= 0.2 and 10 dB squeezing (Dotted curve). For contrast, 
we also show the Heisenberg limit in a unit circle and the 
ideal conditional quantum state in solid ellipse. 



The corresponding covariance matrix V add 
added verification noise is given by 



V 



add 



l 



i-n 



KHpHxl -AC F h/2 



for the 



(98) 



A more summarizing measure of the verification accuracy 
is the uncertainty product of the added noise ellipse with 
respect to the Heisenberg limit, namely, 



U 



add 



|VdetV add = ^ 
n 1 — 7] 



(99) 



In the ideal case with r\ = 0, this simply recovers 
the order-of-magnitude estimate given in the subsection 
IIII Dl In Fig. HI we show the uncertainty ellipse for the 
added noise in the case of ( x = O = 0.2, readout loss 
r] = 1% and with (Green dotted curve) or without (red 
long-dashed curve) 10 dB input squeezing. In compari- 
son, we also plot the Heisenberg limit (unit circle) and the 
conditional state obtained through an ideally noiseless 
state preparation (blue solid ellipse). As figure shows, 
the least challenging scenario already begins to charac- 
terize the conditional quantum state down to the Heisen- 
berg Uncertainty. In this two cases, we have A = 1.48 
and 0.62 respectively, leading to 

C/ add = 0.30 (vacuum) , 0.12 (10 dB squeezing). (100) 



VI. VERIFICATION OF MACROSCOPIC 
QUANTUM ENTANGLEMENT 

In this section, we will apply our protocol to verify 
macroscopic entanglement between test masses in future 
GW detectors, which was proposed in Refs. [HI, [Hj]. In 
the experiment as shown schematically in Fig. 1 1 01 mea- 
surements at the bright and dark port of the interferome- 
ter continuously collapse the quantum state of the corre- 
sponding common and differential modes of the test-mass 



c E i N 
X = X + X 



x d = x E — x N 




FIG. 10: (Color online) A schematic plot of advanced interfer- 
ometric GW detectors for macroscopic entanglement between 
test masses as a test for gravity decoherence. For simplic- 
ity, we have not shown the setup at the bright port, which is 
identical to the dark port. 



motion. This creates two highly squeezed Gaussian state 
in both modes. Since the common and differential modes 
are linear combinations of the center of mass motion of 
test masses in the two arms, namely x c = x E + x N and 
x = x E — i N , this will naturally generate quantum en- 
tanglement between the two test masses, which is similar 
to creating entanglement by mixing two optical squeezed 
states at the beam splitter (5(| [H?} • 



A. Entanglement survival time 

To quantify the entanglement strength, we follow Refs. 
[34ll35l] by evaluating the entanglement monotone — log- 
arithmic negativity defined in Refs. [EH [59j]. It can 
be derived from the covariance matrix for Gaussian- 
continuous-variable system considered here. The bipar- 
tite covariances among (x E ,p E , i N ,p N ) form the follow- 
ing covariance matrix: 



V = 



Vee Ven 
Vne Vnn 



(101) 



where 

Vee = Vnn = 

Vne = V E n = 



(Vxx + *&)/4 (V x c p + V x d p )/2 
(V c + V d )/2 (V c + V d ) 

v y xp 1 r xp II \* pp 1 * pp J 

ft^xx ^xx) I ^ i.^xp ^xjj) 

(V c - V d )/2 (V c - V d ) 

\ y xp y xp } I \ v pp ¥ pp ) 



(102) 



(103) 
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FIG. 11: (Color online) Logarithmic negative Efj as a func- 
tion of the evolution duration te, which indicates how long 
the entanglement survives. The solid curve corresponds to 
the case where Q.f/2-k = 20 Hz and the dashed curve for 
Q.f/2tt = 10 Hz. To maximize the entanglement, the com- 
mon mode is 10 dB phase squeezed at t > te and t < while 
the differential mode is 10 dB amplitude squeezed at t < 
and switching to 10 dB phase squeezed at t > te- 

The logarithmic negativity Ej^ can then be written as 
E N = max[0, - log 2 2a-/h], (104) 



where er_ = ^(£ - \/£ 2 - 4det V)/2 and S = 

detVNN + detVEE — 2detVNE- In contrast to Refs. 
[34l [35| , now the covariance matrix V corresponds to the 
total covariance matrix V tot after the entire preparation- 
cvolution-verification process. For Gaussian quantum 
states, we have [cf. Eqs. [[TO]), (ggj and ([98])] 

V tot = V{t e ) + V add . (105) 

B. Entanglement Survival as a Test of Gravity 
Decoherece 

When te increases, the thermal decoherence will in- 
crease the uncertainty [cf. Eqs. (|48)l and (|105[) ] and 
eventually the entanglement vanishes, which indicates 
how long the quantum entanglement can survive. Sur- 
vival of such quantum entanglement can help us to un- 
derstand whether there is any additional decoherence ef- 
fect, such as Gravity decoherence suggested by Diosi and 
Penrose [49l. I50I] . According to their models, quantum su- 
perpositions vanish within a timescale of H/Eq. Here Eq 
can be (a) self-energy of the mass-distribution-difference, 
namely 

4 Q) = fdxdy G[p(x) - p'(x)][p(y) - p'(y))/r (106) 

with p denoting the mass density distribution and r = 
|x — y|; Alternatively, it can be (b) spread of mutual 
gravitational energy among components of the quantum 
superposition, namely 

4 h) = f dxdyGp(x)p'(y)5r/r 3 / 2 . (107) 



with 5r denoting the uncertainty in location. For the 
prepared test-mass quantum states with width of 5x q , 
we have 

t^ « n q /(G P ) , r<?> « a 1 /^^i/ 2/(Gm 3 /2) . (108) 

where L is the distance between two test masses. Plug- 
ging the typical values for LIGO mirrors with p = 
2.2 g/cm 3 , the separation between two input test masses, 
10m and to = 10 kg, we have 

t£ = 4.3 x 10 9 s, t£ = 1.2 X 10~ 5 s. (109) 

It is therefore quite implausible to test model (a); while 
for model (b), f^r^ is less 0.01 with tt q /2ir = 100 Hz. 
In Fig. 1111 we show the entanglement survival as a func- 
tion of evolution duration. As we can see, the model 
(b) of gravity decoherence can easily be tested, for the 
entanglement can survive for several times of the mea- 
surement timescale T q , which is much longer than the 

predicted Tq . 

VII. CONCLUSIONS 

We have investigated in great details of a followup ver- 
ification stage after the state preparation and evolution. 
We have showed the necessity of a sub-Hciscnberg ver- 
ification accuracy in probing the prepared conditional 
quantum state, and how to achieve it with an optimal 
time-domain homodyne detection. Including this essen- 
tial building block — a sub-Hcisenberg verification, we 
are able to outline a complete procedure of a three-staged 
experiment for testing macroscopic quantum mechanics. 
In particular, we have been focusing on the relevant free- 
mass regime and have applied the techniques to discuss 
MQM experiments with future GW detectors. However, 
the system dynamics that have been considered describe 
general cases with a high-Q mechanical oscillator cou- 
pled to coherent optical fields. To this respect, we note 
that in our results for Markovian systems only depend on 
the ratio between various noises and the SQL, and there- 
fore carries over directly to systems with other scales. In 
addition, the Markovian assumption applies more accu- 
rately to smaller-scale systems which operate in higher 
frequencies. 
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Appendix A: Necessity of a sub-Heisenberg accuracy 
for revealing non-classicality 

As we have mentioned in the introduction part, a sub- 
Heisenberg accuracy is a necessary condition to probe the 
non-classicality if the Wigner function of the prepared 
quantum state has some negative regions, which do not 
have any classical counterpart. 

To prove the necessity, we use the relation between Q 
function and Wigner function as pointed out by Khalili 
[60| . Given density matrix p, th e Q function in the co- 
herent state basis \a) is equal to [3lT loll lo2| 



1 



(a\p\a), 



(Al) 



which is always positive defined. It is the Fourier trans- 
form of the following characteristic function: 



J(/3,/T) = Tr[, 



a i{3a^ 



Pi 



(A2) 



Here & is the annihilation operator and is related to 
the normalized oscillator position x/5x q and momentum 
p/Spg [cf. Eq. (j2"6]) ] by the standard relation 



a = [(x/Sx q ) + i{p/Sp q )}/2. 



(A3) 



If we introduce the real and imaginary parts of /?, namely, 
/? = j3 r + i/3\, characteristic function J can be rewritten 
as 

J(p T , A) = e-&+fiV a Tr[e if, *W s *,)+iPi<fi/Sp t ) fa ( A 4) 

where we have used the fact that e A e B = e A+B e^- A ' B ^ 2 , 
as [A, B] commutes with A and B. Inside the bracket of 
Eq. (|A4|) , it is the characteristic function for the Wigner 
function W(x,p), and thus 



J(/3r,A) 



1 



dx'dp'e 



-(/3 r 2 +/3, 2 )/2 



(2tt) 2 

e -»0r (x'/Sx q )-jft (p'/8 Pq )\\rf x '^p'\ 



(A5) 



Integrating over $ and /?;, the resulting Q function is 
given by 



Q{x,p) = — I dx'dp'e 



*') 2 i (p-p'v 



W(x',p'). 



(A6) 

This will be the same as Eq. |J7|), if we identify 
W rccon (x,p) with Q(x,p) and 



8x1 





5pI 



(A7) 



which is a Heiscnbcrg-limited error. Since squeezing and 
a rotation of x and p axes will not change the positiv- 
ity of the Q function, Eq. (|A6|) basically dictates that 
the reconstructed Wigner function will always be pos- 
itive if a Heisenberg-limited error is introduced during 
the verification stage. Therefore, only if a sub-Heisenberg 
accuracy is achieved will we be able to reveal the non- 
classicality of the prepared quantum state. 



Appendix B: Wiener-Hopf method for solving 
integral equations 

In this appendix, we will introduce the mathematical 
method invented by N. Wiener and E. Hopf for solving 
special type of integral equations. For more details, one 
can refer to a comprehensive presentation of this method 
and its applications by B. Noble [&3j | . Here we will fo- 
cus on integral equations that can be brought into the 
following form as encountered in obtaining the optimal 
verification scheme: 



dt'C(t,t')g(t') = h(t) , t>0. (Bl) 



with 



C(t,t') = A(t-t')+ 

»min.[t,t'] 



E 



dt"B* a (t-t")B a (t' -t"), (B2) 



where a = 1, 2, . . . and B a (t) = if t < 0. 

Assuming that solution to g(t) be a square-integrable 
function in C 2 (— oo,oo), one can split it into causal and 
anticausal parts as: 



g{t) = g+{t) + g-{t) , 

where g~(t) is causal part 



9-(t) 



0, t > 



j(t), t < 

and <7+(t) is the anticausal part of g(t) 

'g{t), t>0 



9+(t) = 



0, t£0. 



(B3) 



(B4) 



(B5) 



This definition enables us to expand the limits of inte- 
gration in (|B"H and ([B2]) to -co < (t,t' ,t") < oo: 



d£ C(t,t')g + (t') = h(t) , t>0, (B6) 

J — oo 

where 

C(t,t')=A(t-t')+ 

/+oo 
dt"[B* a!+ (t-t")B a , + (t'-t")] (+!tn , (B7) 
OO 
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index (+, t") stands for taking causal part of a multidi- 
mensional function in the argument t" . 

Let us first exercise the method in a simple special case 
when B a (t) = 0, Va, this gives a conventional Wiener- 
Hopf integral equation 



+00 



dt' A(t-t')g(t') = h(t) , t>0, 







which can be rewritten as 

r- + OQ 



dt' A(t - t')g+{t') - h(t) 



(B8) 



(B9) 



(+,*) 



Applying Fourier transform in t and the convolution the- 
orem, one gets: 



+00 



dtl 
2^ 



A(tt)g+(Sl) - h(fl) 



(BIO) 



The spectrum of causal (anticausal) function is simply 

P 00 

h(-m= dtg +( ^(t)e iM . (Bll) 



However, this evident relation is not operational for us, 
as it provides no intuition on how to directly get g±(Cl) 
given <?(il) in disposal. The surprisingly simple answer 
gives complex analysis. Without loss of generality, we can 
assume that g(t) asymptotically goes to zero at infinity 
as: Vi : \g(t)\ < e" 70 '*' where 70 is some arbitrary positive 
number, that guarantees regularity of g(Cl) at —00 < CI < 
00. In terms of analytic continuation g(s) of g(Cl) to the 
complex plane s = CI + ij, the above assumption means 
that all the poles of g(s) are located outside its band of 
analyticity —70 < Im(s) < 70. Thus, the partition into 
causal and anticausal parts for g(s) is now evident: 



ffO) = 9+(s)+g-(s) 



(B12) 



where o+(s)(g_(s)) stands for function equal to g(s) for 
7 > 7o(< — 7o) an d analytic in the half plane above (be- 
low) the line 7 = 7o(— 70) 3 - According to properties 
of analytic continuation, this decomposition is unique 
and completely determined by values of g(Cl) on the 
real axis. Moreover, as a Fourier transform of valid C 2 - 
function, it has to approach zero when \s\ — > 00. For 
more general cases, this requirement could be relaxed to 
demand that 00 should be a regular point of g(s) so that 
lim g(s) = const. This allows to include ^-function and 

|s|— >oo 

other integrable distributions into consideration, though 
makes us to add the constant 3(00) to formula (|B12[) as 



3 Functions and g~(s) are, in essence, Laplace transforms of 

g(t) for positive and negative time respectively with only substi- 



tution of variable s 



additional term. For example, for g(t) = e "1*1, a > 
one has the following fourier transform: 



9(s) 



2a 



2a 



+ s 2 (s + ia)(s — ia) 



(B13) 



that has one pole s+ = — ia in the lower half complex 
plane (LHP) and one s_ = +ia in the upper half complex 
plane (UHP). To split /(f2) in accordance with (|B12[) one 
can use well known formula: 



S±0) 



Res[g(s), s±, k 

, x (s-s± tk ) a " 
i s ±,fc} 



(B14) 



where summation goes over all poles {s + .fe} (with a k is 
the order of pole s+,fe) of g(s) that belong to the LHP 
for g+(s) and over all poles {s-^} of g(s) that belong to 
the UHP for g~(s) otherwise, and Res[g(s),s] stands for 
residue of g(s) at pole s. For our example function this 
formula gives: 



~9-{s) 



(B15) 



Using the residue theorem, one can easily show that: 



g+{t) = e- at , for t> 
9-(t) = e° 



(B16) 
(B17) 



Coming back to the equation (|B10[) . assume that func- 
tion A(fl) can be factorized in the following way: 



(B18) 



where a + (_)(0) is a function analytic in the UHP (LHP) 
with its inverse, i.e., both its poles and zeroes are located 
in the LHP (UHP). One gets the following equation: 



a_(0)a+(0)ff + (0) -h(tt) =0 



(B19) 



To solve this equation, one realizes the following fact: for 
any function /, [/(f2)]_|_ = means that / has no poles in 
the LHP. Multiplication of / by any function g- which 
also has no poles in the LHP will evidently not change 
the equality, namely, [<7_(J7)/(f2)]+ = 0. Multiplying Eq. 
(|B19|) by l/a_(f2), the solution reads 



h(n) 



(B20) 



Performing inverse Fourier transform of g + (£l), the time- 
domain solution g+(t) can be obtained. 

Now we are ready to solve Eq. (|B6|) with the general 
kernel in Eq. (|B7[) . Performing similar manipulations, 
one obtain the following equation for g + (CI) in the Fourier 
domain: 



U + J2 B a B*A g + - B a (K~g+)- - h 

\ a / a 



i p. 



= 0, 



(B21) 
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where we have omitted arguments Q of all functions for 
brevity. Since B a is a causal function, B* is anticausal 
and g+ is causal, (B* t g)^ only depends on the value of g 
on the poles of B* a . Performing similar factorization 



A 



(B22) 



with ip+ (V'-) an d (l/V'-) analytic in the UHP 



(LHP), </>+(- 
in the form: 


-0) = 




1 

.9+ = y 

V+ 


" ft " 


i 

+ + T + 



E 



(B23) 



Even though g + also enters the right hand side of the 
above equation, yet (£?*<? + )_ can be written out explic- 
itly as: 



^ g + (n-. k )Re S [B*(n), 



(fi-f2_, fe ) CT * 



Here {f2_fc} are poles of B*(Q) that belong to UHP, and 
therefore g + (Q-^) are just constants that can be ob- 
tained by solving a set of linear algebra equations evalu- 
ating Eq. (|E>23[) at those poles {fi-^}. 



Appendix C: Solving integral equations in Section FVl 

Here we will use the technique introduced in the previ- 
ous section to obtain analytical solutions to the integral 
equations we encountered in the subsections lV Bl and IV CI 

In the coordinate representation, the integral equations 
for ,9i,2 are the following [cf. Eqs. (751) and d))]: 



df 



Cn(M') 


C 12 (M')" 




'91(f)' 







C 2 i(M') 


C 22 (t,f)_ 




.92(f) _ 




h(t) 



(CI) 



where Cjj (i, j = 1,2) are given by Eq. (55)) and we have 
defined ft (f) = Mi/i(0+M2/ 2 (0- Since the optimal 51,2 (t) 
will automatically cut off when t > tf, we can extend the 
integration upper bound Tint to 00. It brings those equa- 
tions into the right shape considered in the Appendix [B] 



In the frequency domain, they can be written as 

[S 11 g 1 ]+ + [Sug 2 ]+ = 0, (C2) 
[52iSi]+ + [5 22 S2]+-f = h, (C3) 
f = (1 - r,)(ny2){e 2 « + 2( F )[G x (G x ~g 2 )-]+. (C4) 

Here S^- are the Fourier transformation of the correlation 
functions C^. Specifically, they are 

T) + (1 - v)e 2g 



Su 
S12 



(1-V)e 2 ^ q 



2(Q + ui m - ij m )(Q - w, 

5*21 = <5i2, 

A 2 , (l-r 1 )(e 2 i + 2( F )tf q 



2[(n + Wm )2+ 7^1(0 - Wm )2+ 7 2 



(C5) 

(C6) 
(C7) 

]' 

(C8) 



Since 5*11 is only a number, the solution to g\ is simply 
~ 9l = -S^[S 12 ~g 2 ]+- (C9) 



In the time-domain, this recovers the result in Eq. (91 
Through a spectral factorization 



we obtain the solution for g 2 . 
1 f 1 



.92 



ft - 5 11 1 52i(5i252)- + r 



Plugging r into the above equation, g2 becomes 



1 



.92 



1 



h + K G x (G* x g 2 )- 



(C10) 



(Cll) 



(C12) 



with k = m 2 f2g£p. A simple inverse Fourier transfor- 
mation gives gi(t) and g 2 (t). The unknown Lagrange 
multipliers can be solved using Eq. (53J). We can then 
derive the covariance matrix V add for the added verifi- 
cation noise with Eq. (55)) . In the free- mass regime, a 
closed form for V add can be obtained as shown explicitly 
in Eq. (98)). 
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